# Load data
data <- readRDS("Data/panel/panel_data.rds")

# Identify participants who appear at least three times
panelists <- data %>%
    group_by(id) %>%
    count() %>%
    subset(n >= 3) %>%
    pull(id)

# Subset to these panelists
data <- data %>%
    subset(id %in% panelists)

# Get global means
global_means <- data %>%
    select(matches("norm_|violence|_therm")) %>%
    skim() %>%
    as.data.frame() %>%
    select(skim_variable, numeric.mean)

## Save
global_means %>% saveRDS("Data/panel/global_means.rds")

# Set bootstrapping parameters
unique_ids <- data$id %>% unique()
sample_size <- unique_ids %>% length()
n_sims <- 2000

# Initialize stash for bootstrap results
boot_results <- data.frame()

# Iterate over simulations
for (i in 1:n_sims) {
    # Generate bootstrap sample
    boot_data <- gen_clustered_boot_sample(data, unique_ids, "id", sample_size)

    # Select needed variables
    boot_data <- boot_data %>%
        select(
            boot_id, survey_time, starts_with("norm_"),
            starts_with("violence_"), ends_with("_therm")
        )

    # Calculate within-individual standard deviations
    within_ind_sd_data <- boot_data %>%
        group_by(boot_id) %>%
        reframe(
            across(
                everything(),
                ~ sd(., na.rm = TRUE)
            )
        )

    # Scale time by month
    within_ind_sd_data$survey_time <- within_ind_sd_data$survey_time / 30

    # Calculate weighted means
    weighted_means <- within_ind_sd_data %>%
        reframe(
            across(
                c(-boot_id, -survey_time),
                ~ weighted.mean(., w = 1 / survey_time, na.rm = TRUE)
            )
        )

    # Append to stash
    boot_results <- rbind(
        boot_results,
        weighted_means %>% mutate(simulation = i)
    )

    # Print progress
    paste0("Simulation ", as.character(i), " complete") %>%
        print()
}

# Export bootstrap results
boot_results %>% saveRDS("Data/panel/bootstrapped_within_ind_SDs.rds")
